library(SDMTools)
wd='/home/jc148322/rob/';setwd(wd)

#read in the asciis & create base.asc
bdall=read.asc('Bdall.asc')
bdwet=read.asc('Bdwet.asc')


#read in the data
pos=read.csv('coord.csv',as.is=TRUE)
threshold=read.csv('maxentResults.csv')
threshold=threshold$Minimum.training.presence.logistic.threshold
frogdata=read.csv('frog_pos.csv', as.is=TRUE)

######a function to determine sig.diff areas - (ImageDiff function edited)
ascDiff=function (tasc, sig.levels = c(0.025, 0.975))
{	
    tasc[which(is.finite(tasc) & tasc <= sig.levels[1])] = 9
    tasc[which(is.finite(tasc) & tasc > sig.levels[1] & tasc <
        sig.levels[2])] = 10
    tasc[which(is.finite(tasc) & tasc <= 1)] = 11
	return(tasc)
}
########
##################################################################################################################
#CREATE THE LAYERS, EDIT THE DATA
#create the sig.diff ascii
bdall[which(bdall<threshold[1])]=0; #set values below threshold to 0
bdwet[which(bdwet<threshold[2])]=0; #set values below threshold to 0
base.asc=bdall; base.asc[which(base.asc>0)]=1
sig.diff=SigDiff(bdall,bdwet)

image.diff=ascDiff(sig.diff,sig.levels=c(0.025,0.975))
image.diff=image.diff*base.asc #gives not sig diff areas a value
image.diff[which(image.diff==0)]=8
######################################
#calculate areas
cs=ClassStat(image.diff,latlon=FALSE) #get the area of sig.diff areas.


######################################
#subset frog occurrence data and cbind bc01 and bc12
frogdata$bc01= extract.data(cbind(frogdata$Easting,frogdata$Northing),bc01)
frogdata$bc12 = extract.data(cbind(frogdata$Easting,frogdata$Northing),bc12)

frogdata=frogdata[,c('Species','bc01','bc12')]
lorica=frogdata[which(frogdata[,1]=='LITLORI'),]; lorica[,2]=lorica[,2]/10
nannotis=frogdata[which(frogdata[,1]=='LITNANN'),]; nannotis[,2]=nannotis[,2]/10

record.time=c('New','New','New','New','New','New','Old','Old','Old','Old','Old','Old','Old','Old','Old')
lorica=cbind(lorica,record.time)

######################################
#edit and simplify the data for biplot
tdata = cbind(pos$bc01,pos$bc12); #
tdata[,1] = round(tdata[,1],1); tdata[,2] = round(tdata[,2]/25)*25 #round bc01 and bc12
tdata = unique(tdata) #find only unique combinations of rounded bc01 and bc12

pos$bdall[which(pos$bdall<threshold[1])]=0; #set values below threshold to 0
pos$bdwet[which(pos$bdwet<threshold[2])]=0;

##################################################################################################################
#CREATE PLOTS AND IMAGES

#make the biplot
png('compare_frogs.png', width=8,height=8,units='cm', res=300,pointsize=5)

        plot(tdata[,1],tdata[,2], xlab='Annual mean temperature',ylab='Annual rainfall',xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        points(tdata[,1],tdata[,2], col='grey', pch=15)

        points(pos$bc01[which(pos$bdall>0)],pos$bc12[which(pos$bdall>0)],col='tan', pch=15)

        points(pos$bc01[which(pos$bdwet>0)],pos$bc12[which(pos$bdwet>0)],col='forestgreen',pch=15)

        points(nannotis[,4],nannotis[,5],col='red',pch=15)
        points(lorica$bc01[which(lorica$record.time=='New')],lorica$bc12[which(lorica$record.time=='New')],col='dodgerblue',pch=15)
		points(lorica$bc01[which(lorica$record.time=='Old')],lorica$bc12[which(lorica$record.time=='Old')],col='blue',pch=15)
		
dev.off()

######################################
#make the image

Colormap=c("grey",colorRampPalette(c("tan","khaki","forestgreen","#003300"))(100))

cols=c('grey','forestgreen','khaki','tan')
pnts=cbind(x=c(429509,446345,446345,429509),y=c(8205059,8205059,8128594,8128594))

png ('image.png',width=3*dim(bdall)[1],height=1*dim(bdall)[2], bg= "white", pointsize=100)
    par(mfrow=c(1,3))

        image(bdall, col=Colormap, ,xlab="", ylab="")
        legend.gradient(pnts,col=Colormap, c(0,1), title='Original')

        image(bdwet, col=Colormap, ,xlab="", ylab="")
        legend.gradient(pnts,col=Colormap, c(0,1), title='New')
		
		image(image.diff, col=cols ,xlab="", ylab="")
		legend.gradient(pnts,col=cols,round(range(sig.diff,na.rm=T),digits=1), title='SIG.DIFF')
dev.off()


